function [nFs, mXY, vShut, bReload] = Rectangle(bAutoLoad)
% GalvoScanControl
% Scan script template: Rectangle
% Description:
% Scan x by y points within a user-defined rectangle. ON and OFF times as
% well as stimulation intensity (in Volts) is user defined. Scanning order
% can be randomized.
%

% Scan scripts should return four variables (in this order):
%  nFs      The sample rate (1/sec between scan points)
%  mXY      A matrix with two columns in the format [X Y], where
%               X  are x-positions
%               Y  are y-positions
%           The unit of this matrix is millimeters, relative to the origin
%  vShut    Vector with same length as mXY where values indicate shutter
%           open/close state (when a shutter is controlled) or the beam
%           intensity (when the output is directed to an intensity
%           modulated laser module, e.g. Coherent CUBE).
%           Note: Values of vShut should not exceed +/- 5V. Values beyond
%           these limits will be fixed at +/- 5V by the software.
%  bReload  Boolean indicating whether the script should be re-loaded each
%           time before executing. Setting this to true (1) will decrease
%           the repetition rate and should only be used for long programs
%           that are not executed at high frequency.

persistent p_sHeight p_sWidth p_sBeamON p_sRand p_sVerticesH p_sVerticesV p_sBeamDwell p_sShutState p_sShutVolt
bReload = 1;

% Default values
if isempty(p_sHeight), p_sHeight = '2'; end
if isempty(p_sWidth), p_sWidth = '2'; end
if isempty(p_sVerticesH), p_sVerticesH = '20'; end
if isempty(p_sVerticesV), p_sVerticesV = '20'; end
if isempty(p_sBeamON), p_sBeamON = '10'; end
if isempty(p_sBeamDwell), p_sBeamDwell = '90'; end
if isempty(p_sRand), p_sRand = '1'; end
if isempty(p_sShutState), p_sShutState = '2'; end
if isempty(p_sShutVolt), p_sShutVolt = '5'; end

% Ask for new parameters (not if the script was re-loaded automatically)
if ~bAutoLoad
    cAnswers = inputdlg({'Height (mm)', 'Width (mm)', ...
        'Dwell points, horizontal', ...
        'Dwell points, vertical', ...
        'Beam ON duration (ms)', ...
        'Beam dwell intervals (ms)', ...
        'Randomize order (0=no, 1=yes)', ...
        'Shutter state (0=off, 1=on, 2=on/off)', ...
        'Intensity (0-5 V)'}, ...
        'Rectangle', ones(1, 9), ...
        {p_sHeight, p_sWidth, p_sVerticesH, p_sVerticesV, p_sBeamON, p_sBeamDwell, p_sRand, p_sShutState, p_sShutVolt});
    if isempty(cAnswers), return, end
    
    % For some reason, closing the input dialog takes forever.
    % Placing this pause here prevents this.
    pause(0.1)
    
    p_sHeight = cAnswers{1};
    p_sWidth = cAnswers{2};
    p_sVerticesH = cAnswers{3};
    p_sVerticesV = cAnswers{4};
    p_sBeamON = cAnswers{5};
    p_sBeamDwell = cAnswers{6};
    p_sRand = cAnswers{7};
    p_sShutState = cAnswers{8};
    p_sShutVolt = cAnswers{9};
end

vSize = [str2num(p_sWidth) str2num(p_sHeight)]; % grid size X-by-Y, mm
nVertices = [str2num(p_sVerticesH) str2num(p_sVerticesV)]; % number of vertices along X and Y directions
nInterval = str2num(p_sBeamDwell) / 1000; % interval between vertex traversals, sec
% Note: nBeamOnDur is a proportion of nInterval, i.e. it cannot be
% larger than nInterval
nBeamOnDur = str2num(p_sBeamON) / 1000; % beam exposure duration, sec
nFs = 1 / nBeamOnDur * 2; % sample rate, 1/sec -> Hz
bRand = str2num(p_sRand); % randomize order, 1=yes, 0=no

% Check that nInterval ./ nBeamOnDur is an integer value
if ~((nInterval ./ nBeamOnDur) == floor(nInterval ./ nBeamOnDur))
    uiwait(warndlg('The dwell interval divided by the ON duration must be an integer value.'))
    return
end

bShutState = str2num(p_sShutState); % Shutter  0= always off, 1=always on, 2=on/off
nShutVolt = str2num(p_sShutVolt);

% Compute grid locations
vX = linspace(0, vSize(1), nVertices(1))';
vY = linspace(0, vSize(2), nVertices(2))';

% Subtract origin (i.e. subtract 1/2 from each side)
vX = vX - (vSize(1)/2);
vY = vY - (vSize(2)/2);

% Assemble all combined vertices
mXY_pos = [];
for i = 1:length(vX)
    mXY_pos = [mXY_pos; repmat(vX(i), length(vY), 1) vY];
end

% Randomize order of vertex presentations
if bRand
    vOrder = randperm(size(mXY_pos, 1));
    mXY_pos = [mXY_pos(vOrder, 1) mXY_pos(vOrder, 2)];
end

% Extend matrix so that temporal resolution is nBeamOnDur
% Before each movement pad mXY with identical prior rows so that
% the beam moves and stays off at the next position for a duration
% of nInterval - nBeamOnDur
nPad = (nInterval ./ nBeamOnDur) * 2; % sec / sec
mXY = [];
vShut = [];
for i = 1:size(mXY_pos)
    %mXY = [mXY; repmat(mXY_pos(i,:), nPad, 1)];
    % Push into last sequence the current position, so that galvos are
    % already positioned by the same the shutter opens
    if ~isempty(mXY)
        mXY(end,:) = mXY_pos(i,:);
    end
    mXY = [mXY; repmat(mXY_pos(i,:), nPad, 1)];
    switch bShutState
        case 0 % shutter always closed
            vShut = [vShut; zeros(nPad-1, 1); 0];
        case 1 % shutter always open
            vShut = [vShut; ones(nPad, 1).*nShutVolt];
        case 2 % shutter toggle between open and closed state
            vShut = [vShut; nShutVolt; nShutVolt; zeros(nPad-2, 1)];
    end
end

% make sure shutter is off at end of sequence
vShut = [vShut; 0; 0; 0];
mXY = [mXY; mXY(end, 1:2); mXY(end, 1:2); mXY(end, 1:2)];

% Plot traces
if 0
    figure
    hold on
    nTd = 1 / nFs;
    vT = 0:nTd:((size(mXY,1)-1)*nTd);
    stairs(vT, mXY(:,1), 'bo-', 'linewidth', 2)
    stairs(vT, mXY(:,2), 'go-', 'linewidth', 2)
    stairs(vT, vShut, 'ro-', 'linewidth', 2)
    legend('X', 'Y', 'Shutter')
    xlabel('Time (sec)')
    grid
end

return
